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Abstract 

The nonequilibrium steady state of a granular fluid, driven by a random 
external force, is demonstrated to exhibit long range correlations, which be- 
have as ~ \jr in three and ~ ln(L/r) in two dimensions. We calculate the 
corresponding structure factors over the whole range of wave numbers, and 
find good agreement with two-dimensional molecular dynamics simulations. 
It is also shown by means of a mode coupling calculation, how the mean field 
values for the steady state temperature and collision frequency, as obtained 
from the Enskog-Boltzmann equation, are renormalized by long wavelength 
hydrodynamic fluctuations. 

I. INTRODUCTION 

Systems of granular particles, like grains of sand or more ideally glass, plastic or metal 
beads, exhibit different flow regimes |l|], depending on the external forcing. A systematic 
experimental study of the rapid or collisional flow regime as compared to the quasi-static, 
slow or frictional regime was first performed by Bagnold J| using an annular shear cell. 
Later a similar but more refined characterization was made in Ref. ||. The possibility of 
coexistence of different flow regimes was observed in an experimental study of flows down 
an inclined chute 

In several more recent experimental studies of the rapid granular flow regime, more mi- 
croscopic properties have been measured. In Ref. M the fluidization behavior of a vertically 
vibrated two-dimensional model granular material has been investigated using high speed 
photography. Patterns at the surface of a vertically vibrated granular layer, analogous to 
Faraday waves in molecular fluids, have been observed in Ref. || and stimulated the interest 
of many theorists 0. An understanding of these patterns through a derivation of e.g. an 
amplitude equation M from the hydrodynamic description of the system, is still lacking, 
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however. In Ref. || the effect of inelastic collisions on the formation of clusters is inves- 
tigated in a system of particles rolling on a smooth surface and driven by a moving wall. 
Finally, Ref. [10] studies the steady state of a vertically shaken granular monolayer, and 
discusses clustering, inelastic collapse and long range order. 

Even rapid flows of model granular materials are poorly understood in general, since 
complicating effects, such as gravity and interactions with boundaries, have to be taken into 
account. If the model granular material consists of spherical grains with a smooth surface, 
collisions between particles can be characterized only by their coefficient of restitution a 
or their inelasticity e = 1 — a 2 . We assume this coefficient of restitution to be a constant, 
independent of the relative velocity between the colliding particles, and refer to the model 
as the inelastic hard sphere (IHS) model. Dissipative collisions complicate the dynamics in 
a nontrivial way: they may cause the system to become unstable, and give rise, for instance, 
to clustering; they create several new intrinsic length scales that might interfere for small 
inelasticity with the system size L, and for large inelasticity with the mean free path l Q . 

By driving an IHS fluid by boundaries or external fields it can reach a steady or oscil- 
latory state. Due to the existence of these new "cooling" lengths, this state is frequently 
inhomogeneous, where the spatial gradients become larger at higher inelasticity. Only for 
small inelasticity the mean free path is well separated from the scale on which the macro- 
scopic fields vary, and a hydrodynamic description [11] through Navier-Stokes or Burnett 
equations is expected to hold. 

In the present paper we investigate the properties of an IHS fluid that is heated uniformly 
so that it reaches a spatially homogeneous steady state. This way of forcing, where a random 
external force accelerates a particle, was proposed by Williams and MacKintosh [I2| for in- 
elastic particles moving on a line. Peng and Ohta [13| performed simulations on a 2D version 
of this model. In two dimensions the model may be considered to describe the dynamics of 
light disks moving on an air table, a system which has been investigated experimentally in 
Ref. |TJ]]. In three dimensions it can be extended to include gravitational and drag forces, 



making it to some extent relevant for gas-fluidized beds |Lo] when hydrodynamic interac- 
tions are unimportant. A similar IHS model with random external accelerations has been 
used by Bizon and Swinney [|T6[ in their computer simulations to test continuum theories 
for vertically vibrated layers of granular material. 

In the present paper we will describe the randomly driven IHS fluid in two and three 
dimensions, and characterize its nonequilibrium steady state (NESS). The single particle 
velocity distribution function in the NESS has been calculated in Ref. [IT]] from the Enskog- 
Boltzmann equation and was shown to be well approximated by a Maxwellian, except for 
an overpopulated tail ~ exp(— Ac 3 / 2 ), where c is the velocity scaled by the thermal velocity 
and A ~ 1/ y/e. Computer simulations of the one-dimensional system of Ref. |12| showed 



the existence of long range spatial correlations in the steady state, which were addressed 



theoretically in in Ref. [|18[]. Here we will give quantitative predictions for long range correla- 
tions 0] in the two- and three-dimensional NESS. Moreover, we extend the mode coupling 
theory of Brito and Ernst |^U| to analyze how long wavelength fluctuations in the NESS 
renormalize the mean field predictions of kinetic theory, and use this theory to calculate the 
renormalized temperature and collision frequency in the NESS. 

To obtain an adequate description of the structures in steady granular flows one does 
not only need the equations of fluid dynamics for the average macroscopic behavior, but also 
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the spatial correlation functions G a b(r), and their Fourier transforms, the structure factors 
S a b(k). Let 5a(r,t) = a(r,t) — (a(r,t)) with (a = n,T,u a ) be the fluctuations of the slowly 
varying fields a(r,t), i.e. the local density, n(r,t), local temperature, T(r,t), and local flow 
velocity u a (r,t) (a = x,y, . . .), around their average values (a(r,t)). Then the objects of 
interest are the correlation functions in the NESS, which are given by the limit, 

G ab (r) = lim i / dr'(5a(r + r',t)5b(r',t)) 

S ab (k) = lim y-^atMMK-M)}- (1) 

Here (. . .) is an average over some initial distribution, Sa(k, t) is the spatial Fourier transform 
of 5a(r,t), and S a b(k) that of G a b(r). Moreover we consider the unequal-time correlation 
functions in the NESS, defined as 

Fab(k,t)= lim V~ 1 {5a(k,t' + t)5b(-k,t')} } (2) 



where F nn (k,t) is the intermediate scattering function pl| . The dynamic structure factor is 
then 

5„ n (fc,^) =ReF nn (fc,z = in + 0), (3) 

where F a b{k,z) is the Laplace transform of F a b{k,t). 

The paper is organized as follows. In section we show how the macroscopic equations 
for granular flow are modified to account for the external driving/heating by the random 
accelerations. Section [TTJ characterizes the noise of external and internal fluctuations, and 
the structure factors and spatial correlation functions are calculated in sections [IV] and |V]. 
The latter section also presents the mode coupling calculations for the temperature and the 
collision frequency in the NESS. Computational details of our molecular dynamics (MD) 
simulations are described in section |VT|, and section |V11| compares our predictions with 
simulations. Some general comments and conclusions are presented in section |V11J| . 



II. MACROSCOPIC EQUATIONS 

Consider a system of inelastic disks or spheres (IHS) (d = 2, 3), driven by a heat source, 
which is described clS db random acceleration, 

d^i Fj i . , , jS 

-jr = — +£M- 4 
at m 

Here F, is the systematic force on particle i — (1,2, ... , N) due to inelastic collisions. If the 
time constant of the heat source is much smaller than the mean free time to between colli- 
sions, then £i(t) can be considered as Gaussian white noise with zero mean and correlation, 



&a(*)£;/9(tf) = ffiiA^ - 0> (5) 

where a,/3 = {x,y, . . .} denote Cartesian components of vectors or tensors. The over-line 
indicates an average over the noise source. It is understood that the ensemble average in Eqs. 



3 



(H) and (0), denoted by the angular brackets, also includes this noise average. To guarantee 
conservation of total momentum the random force has to obey the constraint J2i£i{t) = 0. 
In thermodynamically large systems this constraint gives a correction to (|5]) of 0(1/N), 
which can be neglected. 

The uniformly heated fluid is described by the standard macroscopic equations of fluid 
dynamics, where the temperature equation is supplemented with an additional source term 
m^Q, and a sink term T to account respectively for the heating and the energy loss through 
inelastic collisions: 

d t n + V • {nu) = 

d t u + u ■ Vw = --V • n 

p 

d t T + u- VT = -— (V- J + 11: Vu)-r + mtt, (6) 
an 

where p = mn, u the flow velocity, and \dnT the kinetic energy density in the local rest 
frame of the IHS fluid. The pressure tensor n Q( g = p5 a /3 + 5U a/ 3 contains the local pressure 
p and the dissipative momentum flux 5U a /3, which is proportional to V a up and contains 
the kinematic and longitudinal viscosities v and ui, defined below Eq. ( |Al|) of appendix 
The constitutive relation for the heat flux, J = -kVT, defines the heat conductivity k. 
For small inelasticity the transport coefficients z/, ui, and k are assumed to be given by the 
Enskog theory for a dense gas of elastic hard spheres (EHS) |2"2"fl . 

To lowest order in the spatial inhomogeneities the sink term, representing the energy loss 
through inelastic collisions, is given by |2c 



r = 2 l0 cuT. (7) 

It is proportional to the granular temperature, to the average Enskog collision frequency 
(lo ~ vT, explicitly given in Eq. (|A2|) of the appendix), and to the coefficient of inelasticity 



e = 1 — a 2 = 2g?7o, where a is the coefficient of normal restitution. To explain the term m^Q 
in (^|) we calculate the energy gain of a single particle due to the random force in a small 
time St. This is done by formally integrating (|J) and averaging over the noise source, i.e., 



±m[vHt + St)-vHt)] =±mt*6t, (8) 

where @ has been used. Note that we have defined the granular temperature as twice the 
average random kinetic energy per translational degree of freedom, so that the Boltzmann 
constant &b does not appear in its definition. 

The above equations provide a consistent description for the heated IHS fluid at small 
inelasticity. The energy is not conserved in inelastic collisions, and consequently the tem- 
perature is not a hydrodynamic mode, but a kinetic mode with a relaxation rate oc ^quj. 
Nevertheless, at small inelasticities (70 <C 1), it is consistent to include temperature among 
the slowly changing macroscopic variables, which describe the dynamics of the system on 
time scales t large compared to the mean free time t = and on spatial scales A large 



compared to the mean free path, Iq = voto, where v$ = y2T/m is the thermal velocity. 

At large inelasticities, where e ~ the temperature is a fast kinetic mode, that 

decays on the time scale to, and cannot be included among the slow macroscopic variables. 
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In that case the IHS fluid becomes athermal, and the slow macroscopic fields only involve 
the density and flow field, as is the case in lattice gas cellular automata without energy 
conservation [0,^ . 

Let us consider the decay of temperature in more detail. For a homogeneous state, the 
fluid dynamic equations (|6]) will have as a solution n(r,t) = n, u(r,t) = 0, and T(r,t) = 
T(t), the latter satisfying 



d t T(t) 



(9) 



For long times the system approaches a steady state with a constant temperature, determined 
by m^Q = As uj ~ a/T, we obtain the mean field prediction ( |A2|) , as deduced from 

the Enskog theory, 



Tv, = m 



2/3 



d-i 



(10) 



Further symbols are defined below (XI) in the appendix. To obtain the final approach to 
the NESS, we linearize Eq. (||) around T E in Eq. (|lCi|) . This yields an exponential approach, 
i.e. 



ST(t) = T(t) -T E = 5T(0) exp[-37 u4 



(11) 



In fact 37oCo> can be identified as the decay rate (0) of the long wavelength components of 
the temperature fluctuations, as derived in section [TV] below fl2"5|). The exact time dependent 
solution of Eq. (§) can be obtained implicitly (t as a function of temperature), and reads 



./ 



lT(t) 



2 Tv, 



(12) 



where 



In \x — 1| — - ln(x 2 + x + 1) + v^3 arctan ( 



2x 



V3 



(13) 



and T is the temperature at t — 0. 



III. NOISE CHARACTERISTICS IN THE NESS 

The goal of this paper is to analyze the effects of spatial fluctuations 5a(r, t) with (a = 
n, T, u a ) around the NESS on hydrodynamic space and time scales. As we are dealing with 
fluctuations, we linearize the nonlinear equations around the NESS, with the result (|A1|) 



of appendix 0. Moreover, to extend the average equations to fluctuating equations, valid 
on mesoscopic spatial and temporal scales, we need to calculate the external noise terms 
£ (r,t) and 8 e *(r,t) that contribute to d t u and d t T in (|Al|). These terms originate from 
the random acceleration which enters in the microscopic equations of motion (|J). By 

starting from the microscopic expressions for the momentum and energy density one finds 
that the noise sources are given by the long wavelength components of 
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r(r,*) = iEi<(*)*(r-r i (*)) 

i 

A 77? . a 

These fields are again Gaussian white noise with zero mean and correlations 

~ 4-777/ 

0cx (r)t) ^ex (r , )tO = __ e 3 ( y( p _ r ^ ( y( t _^ j (15) 

as follows from (|5|). 

Next, we argue on the basis of the hydrodynamic equations @ that there exists, close to 
the NESS, a range of hydrodynamic wave numbers k ^> k*, the so called elastic regime, where 
the dynamics of the fluctuations is the same as in a fluid of elastic hard spheres or disks, 
and is driven by internal noise that will be studied next. The validity of the hydrodynamic 
equations (Bj) and (|A2) is restricted to wave numbers k <C 2tt/Iq (to guarantee separation of 



kinetic and hydrodynamic scales), and to k <C 2n/a, where a is the disk or sphere diameter 
(to guarantee that the Euler equations involve strictly local hydrodynamics). So, for the 
existence of an elastic regime in the IHS hydrodynamics, the following constraints must be 
satisfied 

, * (2-K . , 

r«fc«min -,- . (16) 
L to cr J 



Moreover, following McNamara p6 we can distinguish a dissipative regime, klo <C e [typ- 
ically klo < 0(e 2 )}, and a standard regime, kl ^> e [typically klo ~ 0(y/e)], separated by 
a crossover regime, around klo ~ 0(e). In the dissipative regime, dissipation dominates 
compression effects and sound propagation, which are O(kl ), as well as heat conduction, 
which is 0(k 2 ll). In the standard regime, dissipation effects are of the same order as heat 
conduction. As a consequence, the hydrodynamic modes and their propagation velocities are 
those of a fluid of elastic particles, while the corresponding damping rates of heat and sound 
modes still depend on the inelasticity. Only in the elastic regime, klo ^> k*lo ~ O(^), also 
these damping coefficients attain their elastic values. The above argument applies for small 
enough e = 1 — a 2 = 2d r y , where the inequalities ([IB]) are obeyed, and an elastic regime 
exists and is well separated from the dissipative regime. 

In the elastic regime the equations for the macroscopic deviations from the NESS are 
the same as those for a fluid of elastic hard spheres, deviating from thermal equilibrium. To 
describe fluctuating mesoscopic hydrodynamics on these length scales, one can add internal 

noise £ (?",£) and 6 m (r,t), describing the rapid microscopic degrees of freedom. The noise 
strength of the internal fluctuations can be obtained from the fluctuation-dissipation theorem 
27|f23| for the EHS fluid, and is given in Fourier representation by 



2T r — - - i 

V-^Mtfii-k,*) = —k 2 [v(6ap - k a kp) + P^kp] 5(t - 



8kT 



2 



V- l 6™(k, t)e™(-k, f) = ~^k 2 5(t - t'), (17) 
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where u, v\ and k are the transport coefficients for the EHS fluid, and k a is a component 
of the unit vector k = k/k. The effective noise in the heated IHS fluid may therefore be 

described by the sum of external and internal noise, £(k,t) = £ + £ (&,£), with 

a similar expression for 9(k,t). The noise characteristics of £(k,t) and 8(k,t) interpolate 
between two limiting behaviors and the corresponding noise strengths are given by the sum 
of (|I§ and flT?D. 

Having specified the characteristics of the noise sources in the macroscopic equations, 
we conclude this section by summarizing the Langevin-type equations that describe the 
dynamics of the slow fluctuations. To do so, it is convenient to introduce the Fourier modes 
5a(k, t) exp[ik-r] of the linearized hydrodynamic equations (|A1). The mesoscopic equations, 



valid on hydrodynamic space and time scales, i.e. t ^> to and kl <C 1, then take the form 

d t 8a(k, t) = M(fc)5a(fc, t) + f (fe, t), (18) 

where the components of the vector a are labeled with a = {n, T, Z,_L}. Here a = I refers 
to the longitudinal velocity component ui(k,t) = k ■ u(k,t), and a =_L refers to (d — 1) 
transverse components of u(k, t). The matrix M ab with a,b = {n, T, I, _L} is given explicitly 
in (|A3|) , and f a (k, t) is Gaussian white noise with nonvanishing components for a — T, I, J_ 
and correlation function 



V-'Uk, t)f b (-k, f) = C ab (k)S(t - if). (19) 



The noise strength C ab {k) = 5 a bC a b(k) is obtained by taking the Fourier transform of ([151) 
together with (|T7j), and is only nonvanishing for the following diagonal elements 

C ha - 4mT £o + 8nT 2 k 2 _ ATT + 8nT 2 k 2 
dn d 2 n 2 dn d 2 n 2 

GllW = f| + !^ = £ + !^ 

n p p p 

c±lW = «I + ^! = £ + ^!, (20) 

n p p p 
where the NESS condition of Eq. (||), i.e. Y = 2j uT = ml; 2 , has been used. 

IV. STRUCTURE FACTORS 

The equal-time structure factors, introduced in ([I]), obey the equations of motion 

d t S ab (k) = {M ac (k)S cb (k) + M bc (-k)S ac (k)} + C ab (k), (21) 

c 

which follows by formally integrating dl8f) and using (|l^) . The left hand side of (pip vanishes 
since the structure factors do not depend on time in the NESS. The resulting equation can be 
solved by spectral analysis, or numerically. The spectral analysis is summarized in appendix 
[A], where w\ a and V\ a are respectively the a-th component of the right and left eigenvectors 
of the hydrodynamic matrix M, and Z\(k) is the corresponding eigenvalue. 
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Taking then the scalar product of (pip on both sides with left eigenvectors ( |X5| ) of the 
appendix yields 



2>A«(fc)|S«6(fc)M- fe )> = " ZA (fc)+^(fc) • (22) 



Using the completeness relation (|A6| ) and the fact that off-diagonal elements of C a b in ( |20| ) 
vanish, we obtain 

which is the final result for the static structure factors. The time correlation function (fj) in 
the NESS reduces in a similar manner to 

F ab (k, t)=J2 ex P [z x (k)t]w X a(k)v Xc (k)S cb (k). (24) 

Ac 



The above result corresponds to the Landau-Placzek theory |2T| for hydrodynamic correla- 
tions in the NESS. 

To obtain more explicit results we need the explicit forms of eigenvalues and eigenvec- 
tors, which have only been calculated for small k. The eigenvalue equation can be solved 
numerically for any given wave number, and the results are illustrated in Fig. p] for the 
two-dimensional case. Transport coefficients, equation of state p(n,T), and pair correlation 
function at contact are obtained from Enskog's theory for elastic hard spheres or disks. The 
generic features of the spectrum in Fig. [I] are the same as McNamara's case "a = (3 = 0" , 
illustrated in Fig. 6(a) of his study |2B[] on hydrodynamics of granular materials, which 
corresponds to a temperature and density independent heat source. However, neither the 
equation of state, nor the transport coefficients, used in Ref. fl2"5]| , correspond to the heated 
fluid of inelastic hard spheres, used in the present simulations. In the hydrodynamic regime 
(Hq < 1) all eigenvalues are found to be negative for nonvanishing wave numbers (see ap- 
pendix 0). So, all modes are linearly stable. 

With the help of Mathematica the structure factors in the steady state have been calcu- 
lated numerically from Eq. (^T|) with dtS a b(k) = for a given wave number. The resulting 
structure factors are shown by solid lines in Fig. ^, and will be tested against MD simulations 
in section [V1J| . 

Next, we present analytic results for the dissipative regime (kl -C 7o). The eigenvalues 
on the largest spatial scales can be determined as an expansion in powers of k at a fixed 
value of 7o, with the results 

z±(k) = —vk 2 

z ±(k) — TikvD — D s k 2 

z H (k) = -3 lo oo + D H k 2 . (25) 

In later applications the explicit form of the eigenvectors of M(fc) is needed to lowest order 
in k. They read 



S 



w ± (k) = (0,0,0,1) v ± (k) = (0,0,0,1) 

w±(k) = -^=(1, -g(n)T/toi, ±v D /n, 0) « ± (fe) = -j=(l, 0, ±n/v D , 0) 
™ H (fc) = (0, 1, 0, 0) v H {k) = (g(n)T/3n, 1, 0, 0). (26) 

The coefficients g(n), vjj, and D x with A = {_L, H, ±} are calculated in the appendix. In 
the dissipative regime there are two propagating sound modes (A = ±) with a propagation 
speed vd and a damping constant D$- There is a kinetic heat mode (A = H), with a long 
wavelength relaxation rate Zh(0) = — 37oCo>, in agreement with Eq. (|TID . Therefore, on the 
largest spatial scales the temperature deviations have decayed to zero, and temperature 
gradients do not exist; there is no heat conduction. In addition, there are (d — 1) transverse 
velocity or shear modes (A =-L), which are purely diffusive. The corresponding diffusivity, 



D± = v, has the same form as for EHS. In (S7) — ( A10 ) the coefficients are expressed 
explicitly in terms of thermodynamic quantities and transport coefficients. 

In the standard regime, klo ^> 70 (and 70 small) the eigenvalues for shear and sound 
modes are to leading nonvanishing order the same as for the EHS fluid, where the sound 
waves propagate with the adiabatic sound speed vg of the elastic fluid, which is larger than 



the propagation speed vjj in (A7). The damping of the sound and heat modes, on the other 



hand, are larger than in the elastic fluid due to the inelastic collisions. In the elastic regime, 



defined in section pj, where klo ^ y/lo, all transport coefficients are equal to their EHS 
values. 

In summary, the eigenvalue spectrum z\(k) for the uniformly heated IHS fluid is quite 
different from that of the freely evolving IHS fluid, linearized around the homogeneous 



cooling state [p6 |. In the free case, all shear modes (A =_L) and the heat mode (A = H) 
are unstable in the dissipative regime (small k), and propagating modes do not exist for 
klo <C To- Moreover, there exist in this regime a stable diffusive density mode and a kinetic 
temperature mode, which combine into two propagating modes for kl ~ 0(70), where 
crossover occurs from the dissipative to the standard regime. In the heated case, however, 
all modes are linearly stable, and the sound modes remain propagating in the dissipative 
regime down to k = 0. 



V. EFFECTS OF LONG RANGE CORRELATIONS 

In this section we study static and dynamic structure factors and corresponding cor- 
relation functions at the largest spatial scales. Moreover, we show by means of a mode 
coupling calculation, how average properties, which were calculated in section |IJ on the ba- 
sis of a mean field theory (i.e. the Enskog-Boltzmann equation) are renormalized by spatial 
fluctuations. 

The static structure factors have been calculated in (^). In the dissipative regime 
(klo ~ 7o), the relevant eigenvalues ( 
dominant singularity at small wave number of the structure factors S a b(k) ~ 0(1/ k 2 ) in ( |2"3"D 
originates from pairs of transverse modes, where 2z±(k) = —2i>k 2 , and from anti-parallel 
sound modes, where z + (k) + z-(k) = —2Dsk 2 . We start with the transverse structure factor, 
where only the shear modes in (pE|) contribute, and deduce from the equations above, 



25) and eigenmodes (Eq) are discussed below 



The 
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S±{k ) = S±±{k )^—, (27) 

where the relation fl20D has been used for klo <C 7o- The structure factors S a b(k) for a, b 
derive their dominant small-fc behavior from two anti-parallel sound modes and we obtain 
with the help of (H), flU), ©, and at small k, 



where the sound damping constant in the dissipative regime D$ has been calculated in ( |A9|) . 
It depends on the inelasticity. The contribution of the internal noise is subdominant in this 
regime. In a similar manner we find 

S ab (k) ~ B ab S u (k) (k->0), (29) 

where all nonvanishing coefficients labeled (ab) = (11, nn, TT, nT) are listed in Table |I[ The 
remaining structure factors are of 0(1) as k — > 0. 

Next, we consider the spatial correlation functions G a b(r), which are the inverse Fourier 
transforms of S a b(k). The small-k behavior of S a b(k), obtained above, enables us to calculate 
the large r behavior of the spatial correlation functions. The calculations are given in 
appendix ||. One finds for the leading large r behavior in £/iree-dimensional systems, 

yonpu J r 

G±(r)^-^-(- + -^-)l (30) 



16np Ku IDs) r 



and in two-dimensional systems, 



G||(r) ^ G±(r) ^_L(I + _l_) k (|), (31) 

valid for r L, where L is the linear dimension of the system. The subleading large 
r corrections to fl3~T|) are constant terms, independent of r. In the calculations given in 
appendix ||, these constants depend on a cutoff wave vector /c min = 2tt/L, used to evaluate 
the divergent k integrals occurring in the Fourier inversion of S a b(k). To calculate their 



precise values the subleading small-fc corrections to (|27| ) and (|28|) are required. 

The long wavelength behavior of the time dependent correlation function F a b(k,t) in 
([24]) can be evaluated in a similar manner. We quote the result in terms of the Laplace 
transform, F a b(k,z), from which the dynamic structure factor (|||) follows. In the dissipative 
regime (kl 7 ) we find to leading order for small wave numbers, 



(k,z) 



S±±(k) 



z + vk 2 
1 z + i\kv D + D s k 2 



In a similar manner we obtain 
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F ab (k,z)~B ab F ll {k,z), (33) 

where all nonvanishing coefficients B^ are listed in Table |. The dynamic structure factor 
(0) then becomes 

s nn (k, a) ~ \ E ^ + ^Lny + Dik±- (34) 

It contains only Brillouin peaks, coming from the sound modes. There is no central Rayleigh 
peak, because the heat mode is not a slow, but a fast kinetic mode in this regime. In the 
elastic regime S nn (k,Q) has the standard Rayleigh and Brillouin lines of the EHS fluid. 

The existence of long range spatial correlations shows that the NESS is quite different 
from a thermal equilibrium state ||19|| . In fact, the spatial fluctuations also modify ("renor- 
malize" ) the mean field predictions for the averages and the particle distribution functions. 
In appendix a mode coupling calculation is presented to estimate the renormalization 
effects on the average energy per particle E = (1/N) J2i{\ mv i) an d average collision fre- 
quency uj in the NESS, and we recall their mean field values, i.e. Ee = |gH~e and we given 
by (|A^), i.e. u>e oc nx{n)\/TE, where T E is given in (|TT]|). 

As it turns out, the fluctuation contributions, SE and 5u, are finite and well-behaved in 
three dimensions, but logarithmically divergent in the system size L in two dimensions, so 
that d = 2 is the upper critical dimension. The mode coupling calculations of appendix |C| 
yield then in two dimensions, for large L, 

Tness - T E + — In 



4vr V l 

wness - + ^ In (jj—j > ( 35 ) 

where Ce and are calculated in appendix 0. The argument of the logarithm •JoL/Iq is 
an estimate for the ratio of the values for the right (k ~ 7oAo) an d left (A; min = 2tt/L) 
boundaries of the dissipative k range, where the small-fc behavior in ( ^7|) to (p9|) is valid. 
The logarithmic correction becomes only appreciable for large systems with a size L, much 
larger that the so called homogeneous cooling length It = lo/lo, which diverges in the elastic 
limit. Then the renormalization corrections for small inelasticity vanish as 5T ~ eln(eL//o) 
and 5oj ~ e 2 ln(eL// ). Here we have used the relations Ce ~ e and C w ~ e 2 for e — > 0, 
as can be deduced from the results in appendix y. A similar mode coupling theory has 
been recently used in Ref. |20| for freely evolving granular fluids to calculate the long time 



decay of the energy, which deviates from Haff 's cooling law due to inhomogeneities in the 
hydrodynamic fields, and good agreement between theory and simulations was found. 

Before concluding this section we compare the theoretical predictions for the structure 
factors, with and without internal noise in Eqs. ( p0|) and (2T), as shown in Fig. |^. Inspection 



of ( p7D and (|28 ) shows that only the external noise determines their dominant small-A; 
behavior. The question then arises what are the effects of internal noise, and is it meaningful 
to include it in the theoretical description? The answer is affirmative, as we will show 
below. The steady state solution of (EH]) without internal noise clearly behaves at small 
wavenumbers as k~ 2 , but the k independent plateau values, shown by the full theory (solid 
lines in Fig. 0), are missing. These plateau values represent a very short distance correlation 
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~ S(r). Calculation of the plateau value for S a p(k) yields (l/n 2 V){J2i v iaVip) = (T/p)5 a p, 
i.e. the self-correlation term (i = j) in the definition ([!]) of S a p(k), or more explicitly in Eq. 
( fiDD below. Then, addition of this plateau value to the numerical solution of ( ]21| ) without 
internal noise yields the structure factors S™1 , out (fc), shown as dotted lines in Fig. ^]. For the 
transverse velocity structure factor the dotted line coincides with the solid line in Fig. as 
can be shown analytically from (^0|) and (pD). Only the structure factor for the longitudinal 



flow field, 5 , M Vlthout (A;) differs appreciably from S\\(k) in the relevant intermediate regime, 
0.2<A;a<0.5. 

In section |V11| the simulation results will be compared with the theoretical predictions 
with and without internal noise. As it turns out, the comparison shows convincingly that 
the theory with/without internal noise agrees/disagrees with the simulations. 



VI. COMPUTATIONAL DETAILS 

In the two subsequent sections we describe molecular dynamics simulations performed 
to verify our theoretical predictions. In this section we present the computational details, 
before testing our theoretical results against simulations in the next section. 

The model studied here has been extensively used in computer simulations for the freely 
evolving case (no forcing) [^8H30|j23[| , as well as for the randomly accelerated case in one |II 



and two dimensions |T3"]. We consider a system of N inelastic hard disks having diameter a 
in a two-dimensional square cell of length L, with periodic boundary conditions. The disks 
interact via inelastic collisions with coefficient of normal restitution a. For a colliding pair 
(i, j) of particles having equal masses the postcollision velocities are: 

v* =Vi - ^(1 + a)(Vij ■ a)a 

v* = vj + -(1 + a)(v i:j ■ &)&, (36) 

where Vy = Vi — Vj, the asterisk denotes velocities after collision, and er is a unit vector 
along the line connecting the centers of particle j and particle i. 

The energy loss in consecutive collisions, which is proportional to e = 1 — a 2 , is compen- 
sated by a periodic (in time) and instantaneous perturbation of all velocities by a random 
amount. After every time step At, the velocity of each particle is modified according to 

Vi^Vi + (pi, l<i<N, (37) 

where the components of the vectors (fi are taken from a random distribution of zero mean 
and variance <po (i n practice a Gaussian or a flat function of finite support). The time 
step At of this "heating" or "kicking" is chosen much smaller than the mean time between 
successive collisions of a tagged particle (typically a factor 10 4 smaller), so that the system 
under scrutiny reduces to that described in section [TI]. The opposite limit, where At is much 
bigger or comparable to the mean free path, was considered in Ref. |31| (in the presence of 
an additional external damping or drag force), and in that limit clustering was observed. 
The relation between ip and £o, the variance of the noise term in Eq. (|]), can be deduced 
from the energy fed into the system by the kicks. This yields straightforwardly, 
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& - £■ (38) 

Between the heating events, the motion of the disks is free, which enables us to implement 
an event-driven molecular dynamics scheme with a linked-list method ||32|| . The CPU time 
however scales like N 2 because the lists in all cells need to be updated after each heating 
event . 

Four parameters determine the state of the system: the inelasticity e = 1 — a 2 , the 
packing fraction <fi = irNa 2 / (4L 2 ), the reduced lower wave number cutoff k min a = 2-na/L, 
and the heating rate £q. The values of iV investigated in this article vary between 10 3 
and 10 4 , and we shall restrict our attention to high packing fractions for which the use of 
linked lists implies the most significant reduction of computer time. We consider the cases 
of moderate inelasticities (0.6 < a < 1) and complete inelasticity (a = 0). For the latter 
case, inelastic collapse occurs, i.e. the collision frequency involving only a small number of 
correlated particles diverges, as observed first by McNamara and Young [BJj in freely evolving 
fluids of inelastic hard disks. Also in our system, for high enough inelasticity the heating 
seems never sufficient to prevent the inelastic collapse. For a < 0.5, the inelastic collapse 
has been avoided by introducing a slight modification of collision rule (j36|) , as proposed in 
PHI : in each collision, the velocities are first computed according to the standard procedure 
(yi, V2) — > v* 2 )\ the relative velocity v* l2 is then rotated by a random angle smaller than 
a maximum value (typically less than a few degrees), keeping the center of mass velocity 
fixed. Note that this modified collision rule does not change the total energy loss of the 
colliding pair, and does not introduce any spurious drag or forcing on the particles. 

The structure factors in the NESS have been computed for wave vectors compatible 
with the periodic boundary conditions, i.e. of the form (2ir/L)(n x ,n y ). We have obtained 
the density-density structure factors 

Snn{k) = -^/^ eX P(~ zfc ' r *i)) = hi J^ 65 ^ -2 *''"*) )' ^ 

and the velocity-velocity structure factor, defined as, 



n 2 S af3 (k) = y(y^v ia Vjp exp(-ifc-ry)). (40) 



Here the averages are taken in the spatially uniform NESS. The fluctuation Sg a in the 
momentum density, Sg Q (k) = J2i w^ia exp(— ik ■ r^), and those in the flow field are related 
as Sg a = p5u a , where p is the average mass density in the steady state. The second rank 
tensor S a p(k) is isotropic and can be split into a longitudinal and transverse part, 

Scp{k) = k a k p S\\{k) + (Sap - fz«kp)S ± (k), (41) 

where k = k/k. From Eq. (0), it appears that the knowledge of S nn requires the compu- 
tation of an 0(N) quantity. We can rewrite the velocity-velocity structure factor so that 
its computation also increases linearly with the number of particles. For example, for the 
longitudinal part of S a p(k) in (41) we have, 
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)• (42) 

In practice, the different S aa (k) have been computed for every k lying in the disk \k\ < 
fcmax = 67t/o", then averaged over shells of thickness k min = 2tt/L, to achieve better accuracy. 
Moreover, the statistics in the NESS has been increased by averaging over time. Note that 
the above procedure, which gives insight into the microscopic to large scale structure of the 
system, does not require the knowledge of the hydrodynamic (coarse-grained) density and 
velocity fields. 



n 2 S\Ak) 



V 



( v i ■ fe) exp(-tk-ri) 



VII. SIMULATION RESULTS 
A. Approach and characterization of the NESS 

Before addressing the question of the large scale structure of the inelastic fluid, we 
investigate the validity of the macroscopic description given in section [TT] and [V]. The former 
section gives the mean field results for the steady state temperature T E in fliCf ) and collision 
frequency u>e in ( |A2| ) , based on the Enskog-Boltzmann equation. The latter section and 
appendix C show how the long range spatial fluctuations renormalize these mean field values 
and lead to estimates in for the corrections 5T = T — Te and 5u = u — using a 
mode coupling calculation. 

When the system is initially prepared in a configuration having a temperature T different 
from the steady state temperature Te, the time dependence predicted by the mean field 
result (|l2l) is in good agreement with the numerical data. This is shown in Fig. ||| for a 
system with small inelasticity. When the initial temperature To is much larger than Te, the 
heating at short times is dominated by the inelastic dissipation, and Eq. ( |12"D becomes Haff 's 
homogeneous cooling law for a freely evolving system 

m = ! ( 43) 

where to — V^eC^o) is the mean free time in the initial state. This can be seen from the 
asymptotic expansion of the function / defined in (|T3|), 



71 3 

f(x) ~ V% for x — > 00. (44) 

2 cc 

These analytic results for short times are confirmed by MD simulations, as shown in Fig. |j. 
Moreover, for initial temperatures T -C Te, Eq. ([12]) predicts at short times a linear increase 
of T, as in a heated fluid of elastic hard spheres. Our simulations confirm this behavior. 

Fig. [| shows that the measured kinetic energy per particle T is larger than the temper- 
ature Te predicted on the basis of mean field theory. This effect is noticeable albeit small 
in the results reported in Figs. |3| and ^, which correspond to the nearly elastic limit. For 
the densities studied here, we observe (see Fig. ||) that the correction 5T is positive and 
decreases with decreasing inelasticity. The positive excess in temperature is already present 
for small inelasticity (see e.g. Fig. f|), and vanishes as e — > 0. Note that the above results 
are at variance with those reported by Peng and Ohta |13|], who find that T/Te does not 
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depend on a. The measured collision frequency lo is also larger than the Enskog estimate 
cue, as shown in Fig. |5|. The excess 5u increases with increasing inelasticity. 

We first observe that the simulation data in Fig. ||, where bothu > u E and T > T E , cannot 
be explained consistently as a possible over- or underestimation of the IHS pair correlation 
function xihs at contact by its value for elastic disks. On the basis of ( |i~0"|) and ( |X2| ) we note 
that an overestimation of xihs would increase the collision frequency (|A2|) , and decrease 
the temperature flTU|). These trends are at variance with the observations. The discrepancy 
between the measured and predicted temperatures and collision frequencies is more likely 
due to large scale fluctuations, at least for not too large inelasticities (a > 0.6). These long 
range spatial fluctuations renormalize the mean field Enskog values for the temperature Te 
and collision frequency uje by amounts 5T and 5u. The theoretical estimates (^), based on 
mode coupling arguments, show the correct trends for the dependence of these corrections on 
the inelasticity, i.e. ST ~ e and 5u ~ e 2 as e — > 0, and give a rough estimate of the magnitude 
of these terms as illustrated in the inset of Fig. [| We also point out that the system size 
considered in Fig. [| is much too small for the asymptotic theory ( pop to be applicable. For 
instance at a = 0.8, we deduce from the data (lo = 1.1a, L = 60a) in the caption of Fig. 
[I that JqL/Io ~ 4.9. Consequently, the leading asymptotic term ln(7 I/// ) — 1.6 does not 
dominate the full mode coupling contribution (C2) in appendix |Cj, where subleading terms 
of 0(1) have been neglected. The corresponding ratios ^qL/Iq for Figs. |6| (a = 0.92), || 
(a = 0.6) and pTT| (a — 0) are respectively 46, 127 and 297. This predicts for the systems 
in Figs. ||, H and [11] respectively Tness/^e — 1-07, 1.41 and 1.77, whereas the simulations 
yield for the observed values T/Te — 1.05, 1.45 and 1.5, which agrees quite well. The 
good agreement at a = is unexpected, as the theory is constructed under conditions that 
apply at small inelasticity. However, the renormalized values for the collision frequency, as 
predicted by the mode coupling theory, are much too small. We find in the above Figs. [], |8| 
and [ll] respectively for the theoretical values c^ness/^e — 1-003, 1.06 and 1.12, whereas the 
simulations yield uj/ue — 1.05, 1.36 and 22.9. 

For large inelasticity (a < 0.5), the temperature T and collision frequency u depend 
on the maximum angle of the random rotations, used to avoid the collapse singularity, 
as explained in section |VI[ In fact, uj diverges at 9 = (inelastic collapse), in agreement 
with the observations of Peng and Ohta [0. Intuitively, one might expect that a larger 
randomization of the postcollision velocities (larger B), would more effectively destroy the 
correlations leading to inelastic collapse, and consequently would decrease the deviations in 
T and uj from the mean field and the mode coupling predictions. This happens indeed for the 
temperature, but not for the collision frequency. At packing fraction = 0.07, the collision 
frequency is indeed monotonically decreasing with increasing G (see Table [TT|). However, at 
cf) = 0.2 it diverges at = 0, decreases to a finite value at = 5°, and increases again to 
its value at = 10° (see Fig. |5|). We have no explanation for this behavior. 



B. Fluctuations in the NESS 

In this section we analyze the effects of inelasticity on the large distance behavior of the 
fluid. We have computed the structure factors in the spatially homogeneous NESS of the 
inelastic hard disk fluid as explained in the previous section, and have focused either on 
values of a close to the elastic limit, where the theoretical description is supposed to apply, 
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or on values close to 0, in order to test how large deviations from the theory might be. Local 
mean field values, like and we, reach their steady state values rapidly. However the time 
scale needed for the structure factors S a b(k) and the contributions of spatial fluctuations, ST 
and 5u, to reach their steady state values are diffusive, and increase as /c~^ n ~ L 2 with system 
size. We have checked in the simulations that the large scale behavior of the structure factors 
was properly equilibrated before accumulating the data used to compute the averages. 
First of all, we have tested the isotropy of tensor (BO) by checking that the average 



Vi ■ k) (vj ■ k±) exp(tk-Vij) V with k ■ k± = 0, (45) 



vanishes for k values compatible with the periodic boundary conditions. 

In Fig. ^|we show the density-density and the relevant components of the velocity- velocity 
structure factors. For elastic hard disks the plateau values of S aa (k) around ka ~ 2 extend all 
the way down to k — 0. The excess correlations in the dissipative regime {k < 70/Z0), which 
for S±(k) extend up to k ~ y^fo/lo, are characteristic of the randomly driven inelastic fluid. 
Fig. ^ shows that for small inelasticities the agreement between simulations and theory 
is quite reasonable. The structure factors diverge at small scales like k~ 2 , in agreement 
with the theoretical predictions (^)-(^). The packing fraction has been chosen fairly high 
(4> = 0.63) but lower than the two-dimensional random close-packing of monodisperse disks 
(f)RCP — 0.82 ||33|1 , and inside the liquid region of the phase diagram for elastic hard disks. 



In addition to the gain in computer time, such a packing fraction leads to a mean free 
path Zo smaller than the particle diameter a (e.g. Zo — 0.095a for <p = 0.63). Therefore, 
the hydrodynamic regime will hold up to typical particle diameters, enlarging the range of 
wave vectors where comparison between simulations and theoretical predictions is feasible. 
Moreover, for dense systems, a marked density structure is to be expected at the molecular 
scale, especially at wavelengths close to a (ka ~ 2ir). Fig. ^ shows that for a close to 1, 
this structure is indistinguishable from the structure for elastic hard disks. Note that Su(k) 
and S±(k), although quite structureless for ka > 1, show a weak and broad peak correlated 
with the maximum of S nn (k). This feature is more pronounced as the inelasticity increases, 
as shown in Figs. |H] and [I3J. As the molecular structure of the fluid has not been taken into 



account in the long wavelength hydrodynamic approach of section [iy], the present theory 
cannot explain the structure in Figs. [Tt] and [IB], and the structure factors predicted by the 
theory reach a plateau in the elastic regime (kl > y/jo), given by 

Snn(k) > TL 2 T \ T 

T 

S a p(k) — S a(3 , (46) 

where \ T is the isothermal compressibility of the elastic hard disk fluid. For a close to 1, 
the above limiting behavior is observed numerically for Su or S±, and for S nn only in the 
limit of small packing fraction, where the molecular structure disappears. 

When the inelasticity is increased, the structure factors exhibit the same k~ 2 behavior 
at large scale, but the theoretical expressions are less accurate (see Fig. H). However, the 
theoretical curves are based on the Enskog estimate Te for the granular temperature T, 
whereas for the system corresponding to Fig. || T ~ 1.45Te- Our mode coupling theory 
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predicts here Tness — 1.41Te. Fig. || displays the comparison between theory and simula- 
tion when the measured granular temperature is taken as an input for the hydrodynamic 
description. It appears that the large scale correlations (for which the present theory has 
been constructed) are well described by the theory, as long as the temperature is corrected 
from the mean field Enskog prediction (|H]) to the measured value T. In the case of S±, the 
amplitude only depends on the shear viscosity. The good agreement of the amplitude when 
the temperature is rescaled, while keeping for the shear viscosity the elastic hard disk value, 
suggests that the dependence of the shear viscosity on the inelasticity could be attributed 
only to the change in temperature. At the molecular scale, Sy(fc) and S±(k) appear to be 
correlated to the density- density structure factor (see Fig. [T0|) , in marked contrast to the 
elastic situation, where a plateau value would be reached. Such an effect is beyond the scope 
of our hydrodynamic approach, and is currently under investigation. 

Surprisingly, in the case of complete inelasticity (a = 0), the theoretical structure factors 
give a reasonable picture of the large wavelengths correlations in the fluid (especially for the 
longitudinal velocity component Su(k), as shown in Fig. [n]). When the theoretical structure 
factors are deduced from the measured granular temperature T ~ 1.5Te and not from Te 
(our mode coupling theory predicts Tness — 1.77Te), the agreement for S±(k) improves, 
but the mismatch for Su(k) increases (see Fig. O). At small scales, the density correlations 
differ significantly from the elastic ones (Fig. 13) and the velocity structure factors exhibit 
the oscillatory behavior already present in Fig. |H], with peak positions locked in on the 
peaks in S nn {k). As can be expected from Figs. [11] and [L2|, the large scale correlations are 
compatible with the expected k~ 2 law (see Fig. |14|) unlike the results of Peng and Ohta who 
report a k~ 1A asymptotic behavior for a = 0. However, a log-log plot such as Fig. [14] does 
not allow an accurate evaluation of the scaling exponents. The k~ 2 law is better inferred 
from the direct comparison with theory (Figs. [T2] ). 

Before concluding this section we compare the structure factors for the longitudinal flow 
fields, obtained from MD simulations with two different theoretical predictions in Fig. 
obtained by including or excluding internal noise. First observe that all parameters in 
Fig. | and Fig. | are identical, as well as units on both axes. The simulation results for 
Si\ (k) in Fig. |6] are in excellent agreement with the theoretical prediction (dashed line), which 
corresponds to the solid line in Fig. |2| (internal plus external noise). The dotted line (without 
internal noise) for SJ lthont (k) in Fig. |2] disagrees with the simulations in the relevant interval 
0.2 < ka < 0.5. 

Hence, inclusion of internal noise extends the validity of the asymptotic theory to inter- 
mediate wave numbers. 



VIII. CONCLUSION 

We have presented a theory for the large scale dynamics of a granular fluid that is driven 
into a nonequilibrium steady state by a random external force. Our description combines 
the macroscopic equations of motion for the hydrodynamic fields, accounting for energy 
dissipation through inelastic collisions and uniform heating, together with the fluctuating 
forces. The long range character of the spatial correlation functions is determined by the 
small wave number divergence ~ k~ 2 of the corresponding structure factors. This k~ 2 behav- 
ior is typical for systems which combine conserving deterministic dynamics (conservation of 
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particle number and momentum in collisions) with nonconserving noise []34|, thus violating 



the fluctuation-dissipation condition, and is generic for rapid granular flows that are driven 
by external noise. We also draw attention to the analogy of our equations of motion for the 
fluctuating fields to the Edwards- Wilkinson model [^5| that was proposed for growth of a 
granular surface. In that case the dynamic variable is a scalar field, namely the height of the 
surface, which obeys a similar equation of motion as any of the (d — 1) components of the 
transverse velocity field, u± a (k,t), in our case. The only difference is that in the Edwards- 
Wilkinson model there is only nonconserving noise, whereas in our case both nonconserving 
(external) and conserving (internal) noise are present. 

We have tested our predictions for the structure factors against molecular dynamics 
simulations and have demonstrated that there is quantitative agreement over the whole 
wave number range, if internal fluctuations are taken into account. In our two-dimensional 
simulations we have found deviations in the steady state temperature and collision frequency 
that grow with the inelasticity and the system size. For not too large inelasticities (a > 0.6), 
we have explained these deviations in terms of mode coupling effects of the long range 
fluctuations. 

The phenomenological mode coupling theory, proposed in |2(| and extended here to the 
driven IHS fluid, starts from the same basic ingredients as in the case of elastic fluids PS" 



where that theory was used to calculate the long time tails of the velocity autocorrelation 
function and other current-current time correlation functions. For the elastic case the mode 



coupling theory can be derived from the ring kinetic theory in the low density limit [37 



which accounts for dynamic correlations built up by sequences of correlated binary collisions, 
leading to nonlocal effects in space and time. Such collision sequences correct the mean field- 
type Boltzmann or Enskog kinetic equations for the errors induced by the breakdown of the 
molecular chaos assumption. The ring kinetic theory for rapid granular flows of IHS has been 



developed in Ref. ||38|| , but has not yet been used to derive the present phenomenological 
mode coupling theory from the more fundamental kinetic theory, valid in the low density 
limit. 

In detailed balance models, such as elastic hard spheres, dynamic correlations created by 
correlated collision sequences lead to long time tails, which imply that transport coefficients 
in two dimensions diverge as InL for large systems |36| , |37|| . 

Non-detailed balance models, such as IHS fluids, generically exhibit long range spatial 



correlations ||19|| . In randomly driven IHS fluids, as studied in this paper, these correlations 
between densities and flow fields at distant points in the fluid behave as 1/r in 3D and 
lnr in 2D, and already modify (renormalize) the mean field Enskog-Boltzmann values for 
steady state properties, such as the temperature and collision frequency. In 2D systems 
these renormalization corrections, 5T and 5uj, exhibit the InL divergence, which in the case 



of detailed balance models appears only in the transport coefficients ||19|| . 



At larger inelasticity (a < 0.6), molecular chaos is also violated due to the presence of 
short range velocity- velocity correlations (see Fig. [13]), which are beyond our mode coupling 
theory. A detailed investigation of the small scale structure, which for large inelasticity 
clearly deviates from an equilibrium structure, will be reported in a subsequent publication. 
It is surprising that our description, which is based on the Enskog theory and neglects any 
dependence of transport coefficients on inelasticity, even at a = predicts the long range 
structure reasonably well, provided that the temperature is not taken as the mean field 
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Enskog value, but set equal to the value measured in the simulations. 
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APPENDIX A: 



In this appendix we derive expressions for the transport coefficients that govern the 
decay of fluctuations in the steady state. Linearization of the macroscopic equations @ 
around the NESS, defined in (§) and (|T0D, gives the deterministic part of the following set 
of equations: 



d t Sn = -nV • u 
d f u = — 



Vp + uW 2 u + {yi - i/)VV • u + £ 



P 

dn dn 



5T + 9. 



(Al) 



The noise terms £{r,t) and 6(r,t) have been discussed in section [TTI[ The pressure p is 
assumed to be that of EHS, p = nT(l + Vt d xna d /2d), where Vt d = 2n d / 2 /Y{d/2) is the 
(i-dimensional solid angle, and x{ n ) t ne equilibrium value of the pair correlation function 
of EHS of diameter a and mass m at contact. The kinematic and longitudinal viscosities 
v and u h as well as the heat conductivity k are also assumed to be approximately equal 
to the corresponding quantities for EHS, as calculated from the Enskog theory ||22|| , where 
pv = 1] and pv\ = 2r](d — l)/d + ( are expressed in shear viscosity n and bulk viscosity (. 
The collisional energy loss in (|7]), T = 2 r )QtoT, is proportional to the collision frequency 



(A2) 



as obtained from the Enskog theory. In two dimensions we use the Verlet-Levesque approx- 
imation x = (1 — 70/16)/(l — 4>) 2 . 

By taking spatial Fourier transforms 5a(k,t) = J dr5a(r,t)exp(—tk ■ r) in (|A1|) , one 
obtains the mesoscopic equation (113) with the hydrodynamic matrix 



M(k) 



( ikn \ 

■j ug(n)T/n S'Jquj + Dxk 2 ik2p/dn 



V 



ikv\jn 




ikpj pT 




vik 2 






vk 2 ) 



(A3) 



It contains the coefficients 



g[n) 



D T = — 




n d%' 
Xdn 



(A4) 



In the body of the paper we need the eigenvalues z\(k) of the asymmetric matrix M a b(k), 
and its right and left eigenvectors, which are obtained from 
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M(k)w x (k) = z x (k)w x (k) 

M T (k)v x (k) = z x (k)v x (k), (A5) 

where M T is the transpose of M. Here A = ± labels the sound modes, A = H the heat mode, 
and A =_L labels (d — 1) degenerate shear or transverse velocity modes. The eigenvectors 
form a complete biorthonormal basis, i.e. 

Y^ V Xa,( k ) W na(k) = (v X \Wfj,) = 5 Xt , 
a 

J2\w x (k))(v x (k)\=I. (A6) 

A 

Moreover the eigenvalue equation, det[z(k)I — M(fc)] = 0, is an even function of k. Con- 
sequently, M(fc) and M(— k) have the same eigenvalues, which are either real or form 
a complex conjugate pair. So we choose z x {k) = z x (—k) = z x {k). The corresponding 
eigenvectors of M(— k) in case of the sound modes are obtained from the transformation 
{w + (— k), k)} <-> {w_(k), v + (k)}. All other eigenvectors are invariant under the 
transformation k — > —k. 

By setting z(k) = in the eigenvalue equation, one can verify that there are no zero 
crossings at any finite wave number. So, all eigenvalues have a definite (here negative) sign. 
Consequently all modes of the heated IHS fluid are linearly stable. There is no clustering 
instability |23|J2^| and no instability in the flow field |2"3 |, as in the freely evolving IHS fluid. 



In the dissipative regime (klo <C 7o) the eigenvalue equation is solved by an expansion 
in powers of k, and one finds to dominant orders the eigenvalues in the form and 
eigenvectors in the form (p6|). The eigenmodes to dominant nonvanishing order in k are 



listed in (26). There are id— 1) transverse velocity or shear modes (A =-L), which are purely 
diffusive with a diffusivity D± = v\ there are two propagating modes (A = ±) with a speed 
of propagation vd, and sound damping constant Ds, and a kinetic heat mode (A = H) with 
a nonvanishing z#(0). 

For later reference we also express these coefficients in thermodynamic quantities and 
transport coefficients. The speed of sound vd in the dissipative regime (kl <C 7 ) is: 

dp\ 2p 



It satisfies the inequality vd < vs, where vs is the adiabatic speed of sound in the standard 
regime (70 kl < 1), 

4 = 4 + mU) = m. (as) 



An) \pT J \dp 
The damping constant of the sound modes is 

and the dispersion relation for the kinetic heat mode contains the positive constant 

D ^ K ^ ( 1 n ^ \ (AlO) 

H dn 97otup \ X^n dnT J 

The ratios v 2 D jT , v 2 s /T, D s /VT, and D H /\/T are independent of temperature. 
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APPENDIX B: 



In this appendix we calculate the tails of the spatial correlation functions G a b(r), which 
is done by Fourier inversion of S a b(k). Consider first the tensor fields G a p(r) and S a p(k) 
in Eq. (0) with Sa a (r, t) = u a (r, t) (a, /3 — x, y, . . .) being the components of the flow field. 
Both tensor fields are isotropic and can be split into longitudinal and transverse components, 
i.e. 



G 



a(3\ 



r a rpG\\(r) + (8 a p - r a r /3 )G±(r) 

dk r- - *• * 

j^-^exp(tk ■ r) [k a kpS\\(k) + (5 a/3 - k a k p )S L (k) 



(Bl) 



where the small k behavior of S±(k) and Sn(k) is given in (P7|) and (28) . By contracting the 
second line above with f a fp we obtain 

dk , . , , .. . : , .. , (B2) 



G\\(r) 



)d exp(zk ■ r) [(fc ■ r-YS^k) + [1 - ((fe • r) 2 ]S ± (k) 



k a kp) yields in a similar manner an 



Contraction of the second line of ( |B 1| ) with (8 a p 
expression for G±(r). 

In three dimensions the k integral can be performed explicitly and yields 

dk exp(tk ■ r) 1 



(2vr) 3 k 2 
dk exp(zfc • r 



k ■ r) 



Airr 
0. 



(2tt) 3 k 2 

The resulting large r behavior of G a p(r) is given in (BD). 

In two dimensions the integral in ( P2|) over the azimuthal angle yields for large r 



(B3) 



dk exp(tk ■ r) 



(2n) 2 
dk exp(tk ■ r 



[1- (Jfe.f) 



1 

2^ 
1 



00 dk 

min k 



Ukr) 



2tt 



In 



+ 0(1) 



(B4) 



(2tt) 2 k 2 

In two dimensions the k integral diverges for k — * 0. It should in fact be restricted to 
k > fcmin = 2tt/L, which is the smallest allowed wave number when periodic boundary 
conditions are imposed. To obtain the last equalities one needs the small z behavior of the 
Bessel functions J u (z) ~ (z/2) u /T(u + 1) fl39| . Combining these results with (|B2|) yields for 
large r 



Gn(r) ~G. W 



1 1 
- + 



In 



, - , , ■ (B5) 
57rp \p 2D S J \r J 

The remaining spatial correlation functions G a b(r) involving a,b = {n, T}, are scalar fields. 
Their large r behavior is given by 



G n ,(r) 



dk 

W) 1 



exp(ik ■ r)S a b(k) 



TB ab j hx{L/r) (d = 2) 
Ir (d = 3), 



8npD s \ 1/2? 



(B6) 



where B a b{k) is given in Table |. The subleading corrections of O(r ) depend on the cutoff 
k min . To evaluate these terms requires the subleading small k behavior of S a b(k) in Eqs. 
to 
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APPENDIX C: 



In this appendix we present a mode coupling calculation to estimate the contributions 
of the long wavelength fluctuations in the NESS to some quantity h. Examples are the 
particle distribution functions, the energy per particle E, and the collision frequency uo. 
The fluctuations are correlated over large distances, as a consequence of sequences over 
dynamically correlated collisions, the so called ring collisions [i0|, as shown in section by 
explicit calculation of their spatial tails. 

To calculate their contributions to average quantities like h, one may solve the ring kinetic 
equations ||38|| , or estimate these quantities from a more phenomenological mode coupling 
approach, as developed in Ref. ||36|| . The basic assumption made there is that the state 



of the system rapidly decays to a state of local equilibrium, described by the fluctuating 
hydrodynamic fields a(r,t) = {n(r, t), u a (r, t), T(r, t)}. 

For the quantity under consideration this can be implemented by representing h = 
{1/V) J drh(r), and approximating h(r) by its value in local equilibrium, i.e. 

/>ness = ^ / dr<^(a(r))>, (CI) 

where the average is taken over the fluctuating hydrodynamic fields a(r) in the NESS. To 
carry out the average over the fluctuations we expand hi(a) in powers of the fluctuations 
5a = a — (a) around the NESS, yielding 

/^ness - fy« a )) + 7^7 / dr(5a{r)5b{r))A ab 
1 f* dk 

= hi{{a)) + -J j—}S ab (k)A ab , (C2) 

where summation convention for repeated indices has been used. Here A ab is the matrix 
of second derivatives <9 2 /i;(a)/<9a<9a at a = (a). The asterisk indicates that k integrals are 
restricted to the long wavelength range, k < Jo/lo, the so called dissipative range, discussed 
below (|T6|). In this range the structure factors have the form S ab (k) ~ E ab /k 2 on account of 
(Pto(H). 

For dimensionality d > 3 the fluctuation contribution, 5h = /iness - ^(( a ))> i s convergent 
at small k, and gives only small well-behaved corrections to /i;((a)). However, for d = 2, 
the k integral diverges logarithmically at small k (where k ^ 27r/L), and the excess, 5h, is 
given by 

5h ^ AA, no* dA: ^ AA ^ fl*L\ ^ 

ill J27T/L k 47T V Iq J 

Consequently, the fluctuation contribution Sh in 2D systems is a singular function of the 
system size, L, that diverges in the thermodynamic limit. 

We first apply the above results to the energy per particle, E = (1/N) J drei (a(r)), where 
e;(a) = \pu 2 + |nT is the energy density in local equilibrium. From e;(a) the expansion 
coefficients corresponding to /i;((a)) and A ab in ( |C2|) can be calculated, to yield in d = 2: 

1 f* dk 

5E ^i^j 7^)2 \pSim + pSuik) + 2S nT (k)) . (C4) 
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Inserting ( p7|) to ( p9| ) then yields 



47m 



k 2D, I n 



where the coefficient _B n r is listed in Table |. 

For the collision frequency the analog of (|C1|) is to = (1/N) J dr(n(r)ui (a(r))) with 
W|(a) oc nx(n)^/T given in (|^2|). This gives in two dimensions: 

I f* dk 

6uJ -7T 7^2 [Snn(k)A nn + 2S nT (k)A nT + S TT (k)A TT } 
2n J {2tc) z 

* ] 0UJET * [B rm A nn + 2B nT A nT + B TT A TT ] In (^) , (C6) 
8nnpD s V to / 

where the coefficients B ab are given in Table |, and A nn = d 2 (nu) / dn 2 , A nT = d 2 (nu)/dndT, 
and A TT = d 2 {nu)/dT 2 at a = (a). 

The same method can be used to calculate other averages, as well as particle distribution 
functions. For instance, for the single particle distribution function the starting point would 
be 

/nessW = (1/V) J dr(f t («|a(r))), (C7) 
and the above procedure can be applied at once. 
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TABLES 



ab 


B a b 


11 


1 


nn 


n 2 /vjj 


TT 


g 2 (n)T 2 /9v 2 D 


nT 


-g(n)nT/3v 2 D 



TABLE L Coefficients B ab in Eq. (||). 



e 


1.5° 


3.5° 


5° 


45° 


90° 


w/^E 


8.2 


7.4 


7 


4.3 


4.3 



TABLE II. Collision frequency (normalised by the Enskog value) as a function of 0, for a 
totally inelastic system (a = 0) with = 1600 particles, and packing fraction <p = 0.07. 
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FIGURE CAPTIONS 



1. Dispersion relations z\(k)/u> versus ka for <p = 0.4, a = 0.9; the solid lines refer to 
the real parts for A =_L, ±, and H respectively. Dashed lines represent the imaginary 
parts of the sound mode relaxation rates (A = ±). Here, l /a ~ 0.34, 7o0"/Z o ~ 0.14, 
and yfjocr/lo — 0.64. 

2. Structure factors, S±(k) and S\\(k), as obtained from the full theory (solid lines) with 
external and internal noise. The dotted line represents SJ lthont (k) without internal 
noise with the plateau value added (see discussion at the end of section |V1I| ). The 
parameters are a = 0.92 and = 0.63. Fig. |6] shows that the Sn simulation data agree 
much better with the solid line than with the dotted line. 

3. Granular temperature as a function of time for a = 0.92, = 0.078 (la/ a ~ 3.5), and 
N = 1600 particles. The simulation result is compared to the analytical expression (|i~2"D 
(dashed curve). The initial condition corresponds to a fluid-like configuration of elastic 
hard disks. Here At = 3 10~ 5 (mL 2 /To) 1 / 2 ~ 1.5 10~ 3 t and ip = 5.77 10~ 2 (To/m) 1 / 2 . 
Te is the temperature expected on the basis of the Enskog theory (see Eq. (|I0"D). For 
the above parameters, there are on average 3.7 collisions per time interval At in the 
NESS, and T E /T = 9.3. 

4. Time dependence of the granular temperature obtained in the simulation for the two- 
dimensional system of Fig. [| with At = 3 10" 4 (mL 2 /T ) 1/2 , <p = 1.7 10- 3 (T /m) 1/2 . 
Comparison is made with Haff's law fl43|) for homogeneous cooling (dashed curve), and 
with the full solution of Eq. (|i~2"D (long-dashed curve). Here, T E /T = 0.018. 

5. Measured excess temperature, 5T = T — Te, and collision frequency, 8u = uj — co>e, 
versus coefficient of restitution a, for L = 60a, iV = 917 (0 = 0.2, l = l.ler), and 
for two different values of the maximal random rotation angle, B = 5° and O = 10°. 
Inset: comparison with the predictions of the mode coupling theory, Eq. d35|). 

6. Structure factors versus wave vector for a = 0.92, (f> = 0.63 (lo/& — 0.095), and 
iV = 10201. The simulation data (symbols) have been averaged over 10 2 succes- 
sive configurations, separated by a time interval of 20 collisions per particle. Su and 
Sp er p = Sj_ are respectively the parallel and perpendicular parts of the velocity- velocity 
structure factor, defined by Eq. (flUP. Comparison is made with the theoretical expres- 
sions (full, dashed and long-dashed curves) deduced from Eq. ( |23| ) (compare also Fig. 
0). There is a dissipative regime for ka < Iqct/Iq ~ 0.40, and an elastic regime for 
ko ^ v /7o<j/7o — 2.1. Here, T/Te ~ uj/ue — 1-05 whereas the mode coupling approach 
of section |V| predicts T NE ss/^e — 1-07 and c^ness/^e — 1-003. 

7. Same as Fig. ||, with k beyond the dissipative regime. The dashed curve corresponds to 
the density-density structure factor of an elastic hard disk (EHD) system of the same 
size and packing fraction (dashed curve). All results are deduced from MD simulations. 

8. Structure factors for a = 0.6, = 0.55 (lo/cr ~ 0.15), and N = 10201. The lines are 
the corresponding theoretical predictions. The measured temperature T/Te — 1.45 
and our mode coupling theory gives Tness/^e — 1.41. 
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9. Same as Fig. |8] where for the theoretical expressions (lines), the temperature has been 
set to the measured kinetic energy per particle, T ~ 1.45Te. 

10. Same as Fig. |8] beyond the dissipative regime. The theoretical structure factors are 
not displayed. 

11. Structure factors for a = 0, <fi = 0.63 and N = 10201. Comparison is made with 
the hydrodynamic theory (lines). Here, T/Te — 1.5, whereas mode coupling gives 
Tness/Te ^ 1.77. 

12. Same as Fig. [11], where the simulation data (symbols) are compared to the theoretical 
predictions (lines) for which the temperature has been set to the measured kinetic 
energy per particle, T ~ 1.5Te. 



13. Structure factors up to the molecular scale, for the same parameters as in Fig. [11]. S nn 
for elastic hard disks and the same packing fraction has also been plotted (crosses). 

14. log 10 \ S — S min \ versus \og l0 (ka), for the same parameters as in Fig. [11]. Here S min 
denotes the lowest value of structure factor S. The full, dashed and long-dashed curves 
refer to S nn , Su and S± respectively. 
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